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ABSTRACT 


The  Autonomous  Underwater  Vehicle  (AUV)  must  be  able  to  operate  in  various 
shallow  water  sea  state  conditions.  In  order  to  have  a  precise  navigation  and  steering 
system,  and  efficiently  place  charges  on  underwater  mines,  the  AUV  must  be  able  to  sense 
and  overcome  hydrodynamic  forces  which  are  caused  by  waves.  This  thesis  establishes  a 
model  of  sea  state  conditions  based  on  spectral  analysis,  and  uses  the  model  to  predict 
future  knowledge  of  the  sea.  This  prediction  is  determined  by  the  random  white  noise 
output  of  a  discrete,  distal  filter.  The  development  of  the  discrete,  digital  filter  is 
described  herein.  The  Pierson-Moskowitz  (P-M)  spectrum  which  models  seaway 
elevations  using  linear  wave  theory  is  used  as  a  target  spectrum  which  the  filter  will  track. 
Cross-correlation  between  the  P-M  target  spectrum  and  digital  filter  have  shown  that  a 
reasonably  accurate  estimate  of  wave  elevations  can  be  predicted  one  full  wave  period  into 
the  future. 
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L  INTRODUCTION 

A.  fflSTORY  ON  MINE  APPROXIMATION  COUNTERMEASURES 

Once  the  Navy’s  primary  mission  shifted  from  “Sea  Control”  (blue  water 
operations)  to  “From  the  Sea”  (shallow  water  operations),  the  increased  requirement  for 
amphibious  warfare  operations  demanded  that  increased  interest  be  placed  in  the  area  of 
mine  countermeasures  -  specifically,  mine  warfare  in  shallow  water.  Because  of  the 
immense  level  of  tension  and  somewhat  frequent  conflicts  recently  arising  in  the  Middle 
East,  the  interest  of  the  United  States  security  has  often  been  channeled  toward  the  mine 
infested  waters  of  the  Arabian  Gulf  Over  the  past  ten  years  or  so,  the  United  States  has 
discovered  and  cleared  a  large  number  of  underwater  mines  in  this  area.  Furthermore, 
there  have  been  some  adverse  developments  occur  as  an  aftermath  of  ships  encountering 
explosive  mines.  To  this  end,  we  must  continuously  remain  attentive  of  this  phenomenal 
threat  to  our  forces  and  continue  ,  to  seek  improvement  toward  measures  of  combating 
this  hostile  and  unwavering  enemy  by  devising  safe  and  efficient  tools  to  prevent  mines 
from  presenting  future  problems,  and  hindering  the  effective  operation  of  Amphibious 
Warfare. 

As  a  measure  to  alleviate  the  danger  of  underwater  mines,  the  Naval  Postgraduate 
School  Center  for  Autonomous  Underwater  Vehicle  (AUV)  Research  through 
Mechanical,  Electrical  Engineering,  and  Computer  Science  Departments  are  continuously 
researching  new  methods  to  overcome  this  ongoing  challenge. 

B.  MINEHUNTING  IN  SHALLOW  WATER  (10  TO  40  FEET) 

With  the  new  mission  (From  the  Sea)  in  place,  today’s  Modem  Navy,  more  than 
ever  before  has  focused  on  means  to  intelligently  and  efficiently  counter  the  challenge 
presented  by  floating  mines  by  employing  countermeasures  in  shallow  water.  Mine 
hunting  in  shallow  water  has  presented  continuous  challenges  to  the  safe  and  effective 
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operation  of  the  United  States  Naval  Forces.  For  maximum  flexibility  the  mine 
countermeasure  efforts  should  be  covert,  cost  effective  and  relatively  quick,  (Richwine, 
1995).  The  countermeasure  methods  currently  used  by  Naval  Forces  include  the 
involvement  of  humankind,  noisy  platforms  such  as  ships  and  helicopters.  These 
aforementioned  methods  require  a  high  risk  of  human  sacrifice,  and  compromises 
covertness  respectively.  As  another  measure,  marine  mammal  systems  and  special  military 
forces  are  capable  of  operating  covertly.  However,  they  are  scarce  resources  that  require 

extensive  training  pipelines,  (Hunt,  1995). 

Owing  to  the  robust  advancement  in  technology,  countermeasures  using 
autonomous  systems  are  currently  being  built  and  tested  at  the  Naval  Postgraduate 
School.  In  previous  studies  at  the  Navy  Postgraduate  school,  an  autonomous  underwater 
vehicle  (Phoenix)  has  been  developed  and  operated  in  a  test  tanks  where  sea  state 
conditions  were  negligible.  In  order  to  advance  the  studies  and  operability  of  the  AUV , 
wave  and  current  forces  must  be  considered  and  incorporated  in  the  test  environment. 

C.  DISTURBANCE  REJECTION 

The  main  source  of  the  dynamic  forces  encountered  by  underwater  vehicles  or 
submerged  vessels  are  wave  and  current  induced.  These  disturbance  forces  arise  fi'om 
buoyant  and  inertial  effects  due  to  ocean  wave  kinematics.  In  the  past,  most  sea  state 
models  were  predicted  using  an  erogodic  process.  To  best  assess  the  performance  of  an 
AUV  in  shallow  water,  it  is  necessary  to  generate  a  dynamic  model  of  the  velocity  and 
acceleration  particles  of  water  to  model  the  wave  forces.  Although  fi-equency 
decomposition  is  a  well  known  method  to  simulate  time  history  variation  of  the  ocean 
surface  due  to  waves,  the  digital  representation  of  the  sea  surface  elevation  and  pressure 
profile  process  was  thought  to  offer  advantages,  and  is  referred  to  as  a  time  series  of  wave 
history.  Because  there  is  a  strong  uncertainty  of  random  ocean  wave  fields,  in  the  past,  a 
spectral  approach  has  been  best  suited  for  describing  the  intensity  and  frequency 
characteristics  of  ocean  wave  fields. 
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Recent  research  on  simulating  ocean  disturbance  was  conducted  by  Spanos  during 
the  an  Offshore  Structures  Conference,  (Spanos  and  Mignolet,  1986).  His  research 
expanded  on  the  use  of  the  Linear  Prediction  Theory  using  Autoregressive  (AR)  and 
Autoregressive  Moving  Average  (ARMA)  Models.  With  the  former,  standard  method, 
considerable  numerical  difficulties  were  encountered  in  getting  realistic  simulations,  while 
the  latter  methods  experienced  difficulties  accurately  approximating  the  target  (P-M) 
spectrum,  and  an  accurate  approximation  could  only  be  feasible  when  using  a  high  order 
AR  scheme,  (ibid).  Spanos  used  dimensionless  coefficients  in  his  version  of  the  ARMA 
model  in  order  to  utilize  the  filter  for  multiple  wind  conditions  and  sea  states. 

It  was  found  that  the  spectrum  of  the  time  series  generated  by  the  AR  model 
exhibited  considerable  narrowband  fluctuations  with  respect  to  the  P-M  target  spectrum. 
The  fluctuations  are  thought  to  be  associated  with  pole  locations  of  the  transfer  function 
with  respect  to  the  unit  circle.  The  dimensionless  coefficients  used  by  Spanos  resulted  in 
high  fi'equency  models  that  fail  to  match  the  P-M  spectrum.  Additionally,  the  form  of 
dimensionless  coefficients  did  not  properly  illustrate  the  behavior  of  ocean  waves  under 
different  sea  state  conditions,  compared  to  the  nominal  one  chosen  for  the  normalization. 
As  it  turns  out,  the  mathematical  form  of  the  P-M  spectrum  cannot  be  realized  exactly  in  a 
linear  filter  transfer  function. 

D.  OBJECTIVE  OF  THESIS 

Even  though  the  spectral  analysis  of  ocean  wave-induced  motions  have  been 
studied  over  a  period  of  recent  years,  problems  and  limitations  continue  to  exist.  These 
limitations  have  made  it  difficult  to  predict  sea  state  conditions;  and  consequently  have 
made  the  task  of  efficiently  operating  Autonomous  Underwater  Vehicles  challenging. 

For  an  Autonomous  Underwater  Vehicle  (AUV)  to  operate  in  shallow  water  with 
adequate  reliability,  wave  forces  and  their  effects  on  the  AUV  should  be  modeled 
mathematically  with  an  adequate  degree  of  accuracy.  A  major  component  of  the 
mathematical  model  of  these  loads/forces  experienced  by  underwater  vehicles  involves 
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ocean  wave  kinematics,  (Sarpkaya  and  Isaacson,  1981).  In  fact,  over  a  period  of  years 
the  Pierson-Moskowitz  (P-M)  spectrum  has  dominated  the  statistical  description  of  ocean 
wave  kinematics,  (Pierson  and  Moskowitz,  1964).  Thus,  in  analyzing  the  impact  of  these 
forces  on  underwater  vehicles,  it  is  worthwhile  to  generate  a  time  record  of  values  of  wave 
kinematics  that  is  compatible  with  the  P-M  spectrum.  Traditionallythis  has  been  achieved 
by  relying  on  the  harmonic  superposition  method,  (Hudspeth  and  Borgman,  1979, 
Shinozuka,  1972,  Shinozuka  and  Wai,  1979),  which  can  be  cumbersome  and  costly  in 
terms  of  computation  time.  Recently,  with  emphasis  placed  on  improving  the 
computational  time,  Spanos  developed  a  procedure  based  on  digital  filtering  techniques. 
His  method  involved  the  design  of  an  autoregressive  filter  (AR)  and  autoregressive 
moving  average  (ARMA)  filter  that  consisted  a  transfer  function  whose  squared  modulus 
was  a  close  of  the  target  P-M  spectrum,  (Spanos  and  Hansen  1981,  Spanos  1983).  His 
digital  filtering  technique  demonstrated  promise  in  computational  efficiency,  nevertheless, 
it  proved  to  be  a  delicate  task  for  spectral  matching,  (Spanos  and  Mignolet,  1986). 

In  this  thesis,  seaway  dynamics  will  be  modeled  as  a  linear  system  with  a  random 
white  noise  input,  where  the  linear  system  transfer  function  is  approximately  mapped  to 
the  P-M  spectrum.  Further,  this  work  will  use  a  digital  filtering  technique  from  two 
transfer  functions  to  approximate  the  P-M  spectrum.  Initially,  the  effectiveness  of 
Spanos’  AR  and  ARMA  methods  will  be  determined  by  testing  the  stability  through 
eigenvalue  analysis.  These  eigenvalues  are  to  be  modified  as  necessary  to  achieve  required 
stability.  This  model  will  then  be  compared  to  the  output  of  a  digital  filter  generated  from 
the  transfer  function  of  an  eighth-order  system.  Once  the  results  of  the  eight-order 
system  are  compared  to  the  results  of  the  Spanos  model  for  best  mapping  of  the  P-M  and 
pressure  profile  spectra,  at  specified  water  depths,  the  system  (eighth-order  or  Spanos) 
with  the  most  accurate  approximation  of  the  P-M  target  spectrum  will  be  further  studied 
in  a  filter  system  and  used  to  predict  sea  states  and  pressure  responses  at  future  times 
based  on  current  and  previous  measurements.  The  accuracy  will  be  measured  using  cross¬ 
correlation  between  the  target  spectrum  and  the  dynamic  equation/innovator. 
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11.  LINEAR  PREDICTION  THEORY  (LPT) 


A.  HISTORY  AND  FORMULATION  OF  LPT 

Digital  representation  of  sea  surface  elevation  and  pressure  profile  processes  is 
called  a  time  series.  The  probabilistic  description  of  ocean  and  wave  induced  responses 
was  introduced  during  the  1950s,  (Spanos  and  Mignolet,  1986).  Prior  to  the  introduction 
of  the  Linear  Prediction  Theory  (LPT),  the  traditional  techniques  for  generation  of 
simulated  digital  sea  wave  height,  velocity,  and  acceleration  records  were  based  primarily 
upon  the  principle  of  wave  superposition.  For  each  point  in  time  domain,  this  technique 
generally  required  the  computation  of  a  sum  of  harmonic  waves,  (Spanos  and  Hansen, 
1981).  Nevertheless,  because  of  successful  research  efforts  in  the  areas  of  simulation  and 
processing  of  digital  data  in  the  time  domain,  modeling  of  sea  wave  time  history  with  the 
Linear  Prediction  Theory  has  proven  to  be  quite  accurate.  The  LPT  method  provides  an 
efficient  procedure  of  generating  a  collection  of  records  possessing  a  power  spectrum 
which  approximates  a  known  target  spectrum,  (ibid).  In  the  LPT,  the  records  of  sea 
waves  are  approximately  compatible  with  the  Pierson-Moskowitz  (P-M)  target  spectrum, 
and  are  obtained  as  the  output  of  a  recursive  digital  filter  with  a  random  white  noise 
input.  The  problem  of  base  period  repetition  of  wave  records  is  thus  avoided. 
Essentially,  this  technique  simply  realizes  the  approximation  of  the  Pierson-Moskowitz 
(P-M)  spectrum  as  the  output  of  a  digital  filter  to  white  noise  excitation.  More  details  on 
the  development  of  the  P-M  Spectrum  formulation  will  be  addressed  latei. 

When  modeling  the  spectrum,  the  velocities  and  acceleration  particles  for  sea 
wave  simulations  can  be  simulated  utilizing  the  Linear  Prediction  Theory  and  the 
following  relationships  for  the  horizontal  velocity  Pv((o)  and  acceleration  Pa(co)  spectra, 
(Shinozuka,  1972) 

P^{co)  =  a)-\G[(o,z,  H)^ P(o}) 
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and 


pXo3)  =  G>\G{m,z,H)p{(D) 


where  P(co)  is  the  surface  elevation  spectrum,  and  G(o),z,H)  is  a  function  of  (o  depending 
as  well  upon  the  total  water  depth  (H)  and  the  vertical  position  (z)  above  the  sea  floor 
where  the  velocity  and  acceleration  simulations  are  desired.  The  independent  vanable  a 
represents  frequency  in  radians  per  second.  Using  Pv(co)  and  Pa(co)  as  target  spectra,  LPT 
can  be  applied  to  obtain  time  domain  velocity  and  acceleration  records.  This  method  of 
simulation  is  based  on  an  all-pole  recursive  filter,  the  input,  as  stated  above,  of  which  is  a 
random  white  noise  process.  A  one-time  computation  of  the  digital  filter  difference 
equation  coefficients  are  required.  Once  these  coefficients  are  computed  for  a  particular 
target  spectrum,  they  can  be  stored  for  future  usage  as  necessary.  The  generation  of  the 
output  time  domain  record  requires  only  simulation  of  white  noise  digital  samples  and 
application  of  the  difference  equation.  No  computations  of  any  trigonometric  or  other 
functions  are  required  after  the  initial  coefficient  generation. 

The  LPT  was  utilized  to  digitally  simulate  the  time  domain  sea  wave  height 
records.  The  target  P-M  spectrum  was  treated  as  one-sided  extending  over  the  frequency 
range  (0,  cob)-  White  noise  records  were  numerically  simulated  and  used  as  input  to  the 
difference  equation  representing  the  digital  filter.  The  difference  equations  and  methods 
for  determining  the  difference  equation  coefficients  can  be  cited  in  (Spanos  and  Hansen, 
1981).  The  Fast  Fourier  Transform  (FFT)  of  the  digital  time  record  generated  by  the 
difference  equation  was  used  to  calculate  the  response  spectrum  of  the  digital  filter 

output. 

The  following  figures  represent  the  results  obtained  using  LPT.  Figure  1 
illustrates  an  example  record  for  an  upper  frequency  ui*  =  7t/2  rad/sec,  filter  of  order  p  = 
25,  and  wind  velocity  of  40  knots.  Figure  2  illustrates  a  portion  of  a  record  which  is  the 
output  of  a  digital  filter  of  order  50,  upper  frequency  of  (Ot  =  7t/2  rad/sec,  and  wind 
velocity  of  40  knots.  Figure  3  is  the  average  power  spectrum  of  300  records  of  1000  time 
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increments  length  for  an  upper  frequency  cOb  =  n/2  rad/sec  and  wind  velocity  of  40  knots. 
A  sampling  interval  of  T  =  2  seconds  between  pulses  of  the  filter  output  to  white  noise 
input  was  used.  This  procedure  was  repeated  for  simulations  of  lower  order  filters  as 
well.  Noticeably,  a  trend  of  decreased  mapping  between  the  P-M  target  spectrum  and 
power  spectrum  waves  of  the  filter  output  was  obvious.  Figures  4  through  7  illustrate 
records  and  power  spectra  of  outputs  of  filters  of  orders  25  and  40  respectively,  with 
frequency  ojb  =  n/3  rad/sec,  a  sampling  interval  of  T  =  3  seconds,  and  a  wind  velocity  of 
40  knots. 
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NORMALIZED  AUTOCORRELATION 
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Figure  1.  Normalized  autocorrelation,  (Spanos  and  Hansen,  1981). 
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Filter  Order  =  50 


Average  Fast  Fourier  Transform  (FFT)  of  300  LPT-simulated  sea  wave  height 
records,  (Spanos  and  Hansen,  1981). 
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Figure  5.  Power  Spectrum  of  LPT-simulated  sea  wave  height  record, 
(Spanos  and  Hansen,  1981). 
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Figure  6.  Portions  of  a  LPT-simulated  sea  wave  height  record, 
(Spanos  and  Hansen,  1981). 
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B.  MODELING  OF  THE  SYSTEM  AND  REAL  TIME-FORWARD 
PREDICTION 

The  general  block  diagram  model  of  the  proposed  overall  system  including  the 
target  spectrum  transfer  and  the  transfer  function  of  the  dynamic  equation/filter  with  the 
innovator  gains  illustrated  is  depicted  below  in  Figixre  8. 


Figure  8.  Block  diagram  of  system 

where  w(s)  is  the  white  noise  input,  Y{s)  is  the  filtered  estimate  of  the  under  wave 
pressure  response,  Y{s)  is  the  measured  pressure,  So(s)  is  the  error  between  Y(s)  and  Y{s)  , 
and  K’s  are  the  system  gains.  The  independent  variable  5  represents  the  Laplace  variable. 

If  the  following  expression  is  stable  then  the  error  (s)  between  the  filter  estimate  and  the 
wave  response  is  a  bounded  value  as  t—>oo. 


Gjs) 

1  +  KG{s) 
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Y(s),  wave  motion  (i.e.  measured  pressure  signal)  is  assumed  to  arise  from  a  random  white 
noise  driven  system  Avith  a  dynamic  transfer  function  where. 


||G(0)lf  =  S{o)) 

Spanos  developed  a  transfer  function  (in  terms  of  z-transformation)  which  is  a  simple 
version  of  the  more  general  higher  order  general  Autoregressive  Moving  Average 
(ARMA)  model. 


m 


k^\ 


This  transfer  function  in  terms  of  z-transformation  notation  was  used  as  the  output  of  the 
digital  filter  which  is  represented  by  the  symbol  Y(t). 

In  Spanos’  dimensionless  model  the  coefficients  ai  and  b;  were  optimally  chosen  to 
minimize  error  between  the  Pierson  Moskowitz  target  spectrum  equation  and  the  transfer 
function  of  the  ARMA  filter.  Spanos  presented  the  coefficients  in  dimensionless  forms  in 
order  to  use  them  in  various  structural  dynamic  applications.  Nevertheless,  when 
attempting  to  model  a  dynamic  equation  for  the  target  P-M  spectrum,  the  transfer  function 
generated  by  the  coefficients  of  the  dimensionless  model  resulted  in  a  narrow  band,  high 
frequency  spectrum  that  failed  to  adequately  resemble  and  map  the  P-M  spectrum.  The 
high  frequency  results  are  illustrated  in  Figure  9.  Additional  problems  were  encountered 
with  normalization  -  it  was  difficult  to  reconstitute  and  stabilize  the  system,  and  the 
transfer  function  had  unstable  zeros,  more  detmls  on  Spanos’  ARMA  in  folloAving 
chapters. 


Mathematical  calculations  determining  the  error  between  the  P-M  target  spectrum 
(wave  and  pressure  responses)  and  the  dynamic  equation/filter  have  been  performed  to 
verify  the  validity  of  the  system’s  block  diagram.  Figure  10  illustrates  a  modified  block 
diagram  of  the  system  which  emphasizes  how  the  error  (e)  is  measured  schematically. 
Assuming  that  the  seaway  can  be  represented  as  the  output  from  a  white  noise  driven 
filter,  generally,  the  error  (e)  goes  to  zero  when  the  wave  and  pressure  responses  G(s) 
equals  or  matches  the  response  of  the  dynamic  equation  S(a}),  or  simply  : 

8— >0;  ||(j(n>)||  = 
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Figure  9.  Time  Record  of  Spanos’  Dimensionless  ARMA  Model 
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Figure  10.  Block  diagram  of  Forward  Prediction  System,  where  5  =  seaway  internal 

dynamics  and  ^PM  ~ 
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Assuming,  for  the  moment  that  all  white  noise  inputs  are  unity. 


e  =  Y-Y  =  -S*M>  +  G(w-Ke) 


therefore. 


8(1  +  GK)  =  (G-S)w 


s  = 


(G-S) 

(1+ga:) 


w 


From  examining  the  numerator  of  the  above  expression,  it  is  obvious  that  when  G=5  the 
numerator  is  zero,  therefore  s=0,  and  if  GK is  stable,  sis  bounded  and  depends  on 

I|G-5||^1G||-W 

since  G  and  5  are  both  positive,  real. 

The  equation  to  follow  is  given  in  terms  of  a  transfer  function,  and  has  been  used  to 
mathematically  model  the  transfer  function  of  the  dynamic  system.  In  the  work  of  Spanos, 
the  equation  is  given  as  the  digital  filter  of  the  followng  transfer  function  represented  in 
terms  of  z-transformation: 


G(z) 


.-Jk 


Jk=l 


In  this  domain,  z  is  the  complex  variable  z  =  and  s  is  the  continuous  time  Laplace 
variable.  For  frequency  inputs,  5  =  jco  (co  is  the  radian  frequency),  in  terms  of  the 
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corresponding  frequency  response  of  the  discrete  time  system,  we  substitute  z  —  and 
the  magnitude  of  the  system  transfer  function  becomes 


-jwkT\\ 


;=0 


-jwkT 


Jk=l 


which  depends  on  both,  a  and  T. 

For  Time  Domain  Representation: 

Take  G  from  the  model  and  convert  it  to  the  followng  equation  in  the  continuous  time 
domain. 


jc(r)  =  Ax{t)  +  Bu(t),  y  =  Cx 

where  0  <  t  <  co,  and  x(t)  is  a  continuous  state  and  u(t)  is  a  continuous  input  of  the 
interval  0<t<  oo,  with  the  initial  values  x(t  =  0),  u(t  =  0)  and  u(t  >  0)  being  a  series  of 
zero  order  hold  data  such  that  at  4  =  kdt,  u(t0  is  a  data  sequence:  Between  the  intervals 
tk<t<  tk+i,  yvhere  u(t)  =  (Figure  11). 


21 


Figure  1 1 .  Represents  Discrete  Time  Domain 
The  solution  for  x(t=ti^  is  given  by  the  following  integral 

fi+i  ^*+1 

j  x{t)dt  =  J  Ax{t)dt  +  J  Bu{t)dt 

<*  h  '* 

giving 

Xk^ri  =  (e‘^)Xk  +  (e‘^ -  l)a^buk  =>  xi^i  =  4‘^k  +  'yuk 

where 

^  =  e“^  and  y- 

This  converts  x=  Ax +  Bu  Xo  the  following  system  in  discrete  sampled  time  domain 
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x(0  =  Ox(r-i)+rw(0 

yit-l)  =  Cxit-l) 

Initially,  Spanos  used  pole  placement  to  design  the  control  matrix,  K  .  However, 
the  Linear  Quadratic  Regulator  (LQR)  proved  to  be  more  convenient  for  determining  the 
values  of  the  gain  matrix,  leading  to  the  estimator  (  «  in  this  case  is  zero  mean,  white 
noise) 

x,„=t>x,+K(p,-p)  ,  p  =  ci, 

where  pk  is  defined  as  the  measured  pressure  from  the  Pierson-Moskowitz  and  pressure 
profile  spectrum. 

Simplifying  the  above  equation  yields; 

f„=(4.-KC)f.+^.;A=Ci, 

The  best  estimate/prediction  forward  in  time  is  then  using  the  above  model  without 
data  correction  obtained  through  the  folloAving  equation  fi-om  the  discrete  time  domain. 
This  forward  prediction  is  always  based  on  the  estimated  quantity  of  x(/)  ; 

p{t+k)  =  ^'‘Cx{t) 
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C.  PROBLEMS  ENCOUNTERED  WHEN  USING  LINEAR  PREDICTION 

THEORY  (LPT) 

This  method  superseded  the  harmonic  superposition  method.  However,  it  proved 
to  be  quite  cumbersome  for  data  processing  machines  requiring  immense  computational 
time.  Problems  were  noted  with  the  fluctuation  of  the  LPT-generated  spectrum  from  the 
target  spectrum.  Figures  3,5  and  7  illustrate  these  fluctuations.  The  LPT  is  not  always 
computational  eflScient  for  simultaneous  simulation  of  sea  wave  velocity  and  acceleration 
records  at  several  times,  and  points  with  any  combination  of  horizontal  and  vertical 
coordinates,  (Spanos  and  Hansen,  1981).  Considering  the  behavior  of  the  power  spectra 
when  mapping  it  with  the  target  spectrum  as  the  order  of  the  filter  p  was  decreased,  it 
was  discovered  that  increasing  the  filter  order  p,  the  local  frequency  breadth  of  the 
spectral  fluctuation  decreases,  (ibid). 

Using  an  approximation  of  the  P-M  spectrum  as  the  output  of  digital  filters  to 
white  noise  excitation  has  been  proven  to  have  numerical  difficulties  when  approximating 
it  by  the  Autoregressive  (AR)  models,  (Spanos,  1983).  However,  when  using  the  AR 
model  as  an  initial  spectral  approximation,  it  proves  to  yield  an  efficient  autoregressive 
moving  average  (ARMA)  model,  (Spanos  and  Mignolet,  1986). 
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HL  APPLICATION  TO  WAVE  MOTION  PREDICTION 


A.  PIERSON-MOSKOWITZ  (P-M)  SPECTRUM 

A  semi-empirical  expression  for  the  frequency  spectrum  of  fully  developed  seas  is 
given  by. 


4 


This  expression  was  proposed  by  Pierson  and  Moskowitz  for  the  height  of  fully 
developed  seas,  (Pierson  and  Moskowitz,  1964).  The  a  and  >9  coefiBcients  are 
nondimensional  parameters  defining  the  spectrum.  The  P-M  spectrum  was  developed  for 
vwnd  velocity  at  a  standard  height  of  19.5  meters  above  the  free  surface.  This  two 
parameter  spectrum  is  sufficiently  general  to  fit  quite  a  few  observations  and  is  consistent 
with  theoretical  predictions  of  the  high  frequency  limit.  Where  most  common  values  for 
the  nondimensional  parameters  are 

a=  8.1x  10'^ 

)9=-0.74 

The  gravity  constant  is 

g=  32.174  fl/sec^ 

The  two  variables  in  the  expression  are 

(o  =  frequency  (rad/sec) 

V  -  wind  speed  (ft/sec) 


The  Pierson  and  Moskowitz  expression  can  be  further  reduced  by  considering  the 
relationship  between  the  significant  wave  height,  Hy}-,  and  wind  speed,  V.  The  following 
expression  represents  the  aforementioned  relationship: 
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=  0.2092 — 
g 


H, 


1/3 


Substituting  this  relation,  along  with  the  specified  nondimensional  parameters  and  gravity 
constant,  into  the  Pierson  and  Moskowitz  spectra  equation  yields  a  simplified  version  of 
the  P-M  expression  which  was  used  to  generate  the  P-M  and  pressure  profiles  for  the 
target  spectrum  (S((o)).  The  simplified  expression  is 


v  8.1x10-^x32.174^ 

S(co)  = - i - X  exp 


Q) 


-0.74 


^32.174^  X  0.2092^^ 


where  H  was  for  (the  mean  water  depth)  shallow  water,  linear  waves  for  this  particular 
model  is  7  meters  (approximately  23  feet).  (Papoulias,  1993) 

Typical  plots  of  varying  significant  wave  heights  (Hs)  versus  fi-equency  and  period 
for  P-M  spectra  are  illustrated  in  Figures  12  and  13  respectively. 


B.  PIERSON-MOSKOWITZ  (P-M)  AND  PRESSURE  PROFILE  TARGET 

SPECTRUM 

In  performing  digital  time  series  analysis,  the  sampling  period  T  and  the  cut-off 
fi'equency  <u,  are  related  through  the  Nyquist  relation  T  =  In  order  to  generate  the 
appropriate  target  spectrum  at  a  specified  operating  depth,  in  this  case  5.2  meters 
(approximately  17  feet)  below  the  surface,  the  pressure  at  this  depth  had  to  be  considered. 
To  this  end,  the  P-M  expression  was  multiplied  by  G((o,z,H)  which  resulted  from  the  linear 
wave  theory.  Plots  consisting  of  varying  significant  wave  heights  (Hs)  versus  fi-equency 
and  period  for  P-M  and  pressure  profile  spectra  are  illustrated  in  Figures  14  and  15 
respectively. 
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Figure  12.  Pierson-Moskowitz  (P-M)  spectra  for  significant  height  (Hs)  of  8,  6,  4,  and  3 

feet  versus  frequency. 


Period-T  (seconds) 


Figure  13.  Pierson-Moskowitz  (P-M)  spectra  for  significant  height  (Hs)  of  8,  6,  4,  and  3 

feet  versus  periods. 
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S(w)  meai 


Figure  14.  Pierson-Moskowitz  (P-M)  and  pressure  profile  spectra  for  significant  height 

(Hs)  of  8,  6,  4,  and  3  feet  versus  frequency. 
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S(w)  mean  square  feet  per  rad/sec) 


Figure  1 5 .  Pierson-Moskowitz  (P-M)  and  pressure  profile  spectra  for  significant  height 

(Hs)  of  8,  6,  4,  and  3  feet  versus  periods. 
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The  Depth  Modifier  for  pressure  variation  expression  is 


cosh(2n(7/  - 

Pwave(t,X,z)  p§Z"H^aP8 

cosh(2n^) 

sin 

2nf---l 
Vt  x)\ 

Revised  Depth  Modifier  for  the  pressure  variation  expression  neglecting  gravity 


G(z,H,X)  = 


cosh(2n(/f 


The  expression  for  generating  P-M  and  pressure  profile  target  spectrum  is 


8.1x10"'  X  32.174" 


0' 


X  exp 


-0.74 


32.174"  X  0.2092" 


V 


X 


cosh(2n(.^  - 


Taking  the  regular  wave  theory,  and  assuming  a  horizontal  sea  bottom  and  a 
fi-ee-surface  of  infinite  horizontal  extent,  a  linear  wave  theory  for  propagating  waves  can 
be  derived.  This  derivation  applies  the  fi'ee-surface  and  sea  bottom  conditions,  (Faltinsen, 
1990),  and  the  Laplace  equation  V"^  =  0.  The  mean  water  depth  will  be  defined  by  the 
constant  H.  The  following  expression  then  gives  the  connection  between  wavelength  X 
and  wave  period  T  for  a  finite  water  depth. 


Substituting  o)  =  27!/Z  in  the  above  expression  and  simplifying  the  nonlinear  quadratic 
equation  yields 


f= 


2n 


=  0 


where  all  values  of  lambda  A  were  found  as  functions  of  frequency  a,  for  this  particular 
depth  condition  see  Table  1 . 

For  this  particular  research,  the  wavelength  and  frequency  connection  equation 
was  solved  for  a  specified  mean  water  depth  (H)  of  7  meters.  The  quadratic  was  solved 
using  the  recurrence  equation  solver  in  conjunction  with  the  floating-point  arithmetic 
solver  from  the  Maple  Software  Application.  The  frequency  range  for  lambda  A  as  a 
function  of  o  is  0.2  <  0)(i)  <  3,  in  increments  of  0.05. 

Thorough  details  regarding  the  derivation  of  linear  wave  theory  for  propagating 
waves  are  beyond  the  scope  of  this  work,  however,  information  can  be  found  in  many 
textbooks  in  fluid  mechanics  (for  instance  Newman,  1977,  chapter  6),  (Faltinsen,  1990). 
The  equations  for  free-surface  and  sea  bottom  conditions,  along  with  the  relations  of 
wavelength  A  and  period  T  can  be  cited  in  (Faltinsen,  1990). 
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171.61 


146.52 


127.62 


112.86 


100.98 


91.21 


83.01 


76.02 


70 


64.73 


60.08 


55.94 


52.23 


48.88 


45.83 


1.05 

43.05 

1.1 

40.49 

1.15 

38.14 

1.2 

35.96 

1.25 

33.93 

1.3 

32.05 

1.35 

30.29 

1.4 

28.64 

1.45 

27.10 

1.5 

25.66 

1.55 

24.30 

1.6 

23.03 

1.65 

21.84 

1.7 

20.71 

1.75 

19.66 

1.8 

18.67 

1.85 

17.75 

1.9 

16.88 

1.95 

16.07 

2.0 

15.30 

14.59 


13.92 


13.29 


12.70 


12.15 


11.63 


11.15 


10.69 


10.26 


9.86 


9.47 


9.11 


8.77 


8.45 


8.15 


7.86 


7.59 


7.33 


7.08 


6.85 


Table  1.  Lambda  (X)  defined  as  fiinction  of  omega  (co)  for  water  depth  of  77  =  7  meters 
(approximately  23  feet). 
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C.  SPANOS’  DIMENSIONLESS  ARMA  FILTER 

Recently,  Spanos  developed  a  digital  filtering  techniques  that  involved  the  design 
of  an  autoregressive  moving  average  (ARMA)  filter  having  a  transfer  function  whose 
squared  modulus  was  assumed  to  be  a  close  approximation  of  the  target  P-M  spectrum. 
The  transfer  function  of  the  AKMA(l,m)  process  y  is  a  discrete  random  process  whose  n* 
sample  can  be  obtained  fi^om  the  previous  one  using  the  following  expression: 

/  m 

k=\  j=0 

where  a  denotes  a  bandlimited  [-a)b,(Ob]  white  noise  process  with  autocorrelation 
function,  (Spanos  and  NDgnolet,  1986). 

Spanos  also  used  the  process  y  as  the  output  of  a  digital  filter  of  a  transfer  function  in  the 
terms  of  z-transform  notation.  The  input  to  this  transfer  is  a  discrete  band  limited  white 

noise  process  o. 


m 

ZV' 

1  + 

ifc=l 

This  transfer  function  was  examined  using  the  dimensionless  coefficients  in  Table 
2.  After  comparing  the  power  spectrum  of  the  transfer  function  and  developing  a  wave 
record,  it  was  discovered  that  this  process  resulted  in  a  very  high  frequency  record  which 
was  far  fi'om  mimicking  the  behavior  of  ocean  waves.  This  behavior  was  cited  earlier  in 
Figure  9.  Consequently,  another  attempt  to  utilize  this  transfer  function  for  matching  the 
P-M  Spectrum  was  made  by  modifying  the  numerator  in  order  to  best  represent  wave 
motions. 
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Manipulating  the  latter  transfer  function  by  taking  the  z  transforms  to  get  a  z-domain 
(ARMA)  model  of  system  dynamics  yields 

zX  =  (1>X  +  yU  =>  (z  -  (1))X  =  yU 
where  ^  e  i?"*"  and  e  i?"*' 

therefore, 

where  y  and  ^  are  m  and  n*-order  polynomial  expressions. 

Spanos  ARMA  model  transfer  function  is  simplified  to  a  more  general  higher  order 
ARMA  model  expression  such  as: 


G(z)=  - 


k^\ 


in  which  bj  and  at  are  related  to  ^'and  (f>. 
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Number 

Coefficient  (b;) 

Coefficient  (a,) 

0 

8.1047x10*^ 

- 

1 

-2.5140x10'^ 

-9.2280x10'* 

2 

1.11 36x1 0'^ 

2.3080x10'^ 

3 

3.7368x10'* 

1.3286 

4 

-3.8057x10'* 

-8.9433x10'* 

5 

-1.2058x10'* 

2.1836x10'* 

6 

2.3579x10'* 

3.4409x10'* 

7 

-8.5629x10'^ 

-1.9221x10'* 

8 

-4.1185x10'^ 

8.6611x10'^ 

Table  2.  Spanos’  P-M  Wave  Elevation  Spectrum  ARMA  Coefficients  (dimensionless). 
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D,  EIGHTH-ORDER  FILTER 

Taking  into  account  Spanos  used  an  eighth-order  transfer  function  to  generate  the 
system  of  dimensionless  coefiBcients,  a  eight-order  filter  was  assumed  to  be  best  suited  for 
generating  a  filter  to  model  the  target  spectrum  in  this  work.  However,  spectral  densities 
or  wave  spectra  for  transfer  functions  of  different  orders,  ranging  fi-om  second  to  eight- 
order  systems  were  evaluated  and  compared  to  the  target  wave  spectrum.  After 
evaluating  the  spectrum  of  all  aforementioned  ordered  systems,  the  spectrum  of  the 
eighth-order  system,  as  expected,  most  closely  reflected  that  of  the  target  wave  specrtum. 
The  comparison  between  the  target  wave  spectrum  and  the  eighth-order  system’s 
spectrum  is  depicted  in  Figure  16. 

The  following  expression  was  used  to  formulate  the  eighth-order  transfer  function  for  the 
dynamic  equation  which  was  used  as  filter  to  forward  predict  seaway  surface  elevation. 


Generally,  the  eighth-order  filter  process  .P  is  a  discrete  random  process  whose 
sample  can  be  obtained  from  the  previous  sample  k  in  the  following  way 

^t+i  ~  A  “  ^^k 

where  the  state  matrix  (p  is  determined  by  the  transfer  function  and  ^  is  a  random  white 
noise  input  with  zero  mean. 
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The  filter  innovator  (gain  matrix  K)  for  this  system  was  generated  by  using  the 
LQR  for  determining  the  stability  of  the  system.  The  fi-equency  was  chosen  according  to 
the  peak  fi'equency  of  the  target  wave  spectrum.  This  value  was  determined  to  be 
approximately  1 .3 1  rad/sec.  It  is  worthwhile  to  note  that  the  peak  fi-equency  for  the  P-M 
and  P-M/pressure  profile  spectra  was  approximately  1.2  rad/sec.  This  frequency  value, 
even  though  the  difference  is  relatively  small  when  comparing  to  that  of  the  P-M  spectrum 
alone,  didn’t  match  the  target  spectrum  as  well  as  the  fi-equency  of  the  P-M  spectrum. 

Tests  conducted  with  the  dynamic  system  using  the  two  different  peak  fi-equencies 
were  conducted.  The  fi-equency  fi-om  the  P-M  spectrum  predicted  sea  condition  out  to  a 
further  time  with  better  accuracy.  The  filter  gains  for  the  LQR  system  were  generated  to 
match  the  target  wave  spectrum.  Noticeable  was  the  small  difference  between  the  filter  of 
the  pressure  and  target  spectrum. 

For  this  work,  significant  wave  height  of  3  feet  was  used  for  both  the  P-M  target 
spectrum  and  pressure  profiles.  Figures  17  through  20  illustrate  the  frequency  and  period 
plots,  depicting  the  significant  wave  height  at  this  particular  depth,  for  each  spectrum. 
Figures  21  and  22  illustrate  the  relationship  between  the  P-M  spectrum  pressure  profiles 
by  normailmng  and  superimposing  them  on  the  same  plot.  Figures  23  through  25 
represent  time  records  of  the  P-M  spectrum,  pressure  profiles,  and  eight-order  filter 
spectrum  respectively.  When  comparing  figure  25  to  figure  9  (Spanos  dimensionless  time 
record),  it  is  obvious  that  they  differ  immensely.  The  a;  and  b;  coefficients  of  the  eighth- 
order  filter  are  listed  in  Table  3. 
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Table  3.  Eighth-Order  Filter  Coefficients. 
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S(w)  mean  square  feet  per  rad/sec 


Figure  17.  Frequency  plot  of  P-M  Wave  Spectrum  for  significant  wave  height  (H  1/3) 

of  3  feet. 
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S(w)  mean  square  feet  per  rad/sec 


Figure  18.  Period  plot  ofP-M  Wave  Spectrum  for  significant  wave  height  (H  1/3)  of  3 

feet. 
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P(w)  mean  square  feet  per  rad/sec 


Frequency-w  (rad/sec) 


Figure  19.  Frequency  plot  of  Pressure  Profile  Spectrum  for  significant  wave  height 

(H1/3)  of  3  feet. 
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P(w)  mean  square  feet  per  rad/sec 


Figure  20.  Period  plot  of  Pressure  Profile  Spectrum  for  significant  wave  height  (H1/3) 

of  3  feet. 
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S(w)  mean  square  feet  per  rad/sec 


Figure  2 1 .  Frequency  plot  P-M  spectrum  superimposed  on  Pressure  Profile  Spectrum 
for  significant  wave  height  (H1/3)  of  3  feet  (normalized  spectra). 
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S(w)  mean  square  feet  per  rad/sec 


Figure  22.  Period  plot  P-M  spectrum  superimposed  on  Pressure  Profile  Spectrum  for 
significant  wave  height  (H1/3)  of  3  feet  (normalized  spectra). 
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Figure  24.  Time  History  of  Pressure  Profiles  reflecting  random  output  for  significant 

wave  height  (Hj/s )  of  3  feet. 
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Wave  Height(ft) 


Figure  25.  Time  History  of  Measured  Pressue  Profile  reflecting  random  output  of 
eighth-order,  digital  filter  for  significant  wave  height  (H1/3  )  of  3  feet. 
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IV.  RESULTS 


To  aid  in  the  discussion  of  results,  it  is  helpful  to  give  the  following, 
Definition: 


p{t,{t  +  t)):  estimate  of  forward  pressure  at  time,  (t+x)  predicted  at  current  time,  (t). 
y(t):  measured  pressure  at  time,  (t). 

y(t  +  x):  measured  pressure  to  be  found  fi-om  seastate  at  time,(t  +  x). 

Prediction  error,  (So)  at  time,  (t  +  x)  based  on  filter  data  up  to  time  (t): 

e[t,  (t  +  t))  =  p[t,  (t  +  x))  -  y(t  +  x) 

Figure  26  illustrates  how  this  error  is  defined  graphically  while  Figures  27  through  38 
graphically  reflect  time  history  plots  of  the  results  obtained  fi^om  the  above  definitions. 

A  measure  of  effectiveness  is  the  forward  prediction  using  the  correlation 
coefficient  between  digital  output  of  the  P-M/pressure  profile  target  spectrum  and  transfer 
function  of  eighth-order  system.  This  measure  of  effectiveness  predicted  wave  and 
pressure  responses  out  to  designated  number  of  time  steps  (number  of  time  steps  for  each 
measure  is  reflected  on  plots  of  the  results).  However,  in  this  work,  results  were  taken  at 
time,  x=5  and  x=30  seconds.  A  fluctuation  was  observed  with  the  cross-correlation 
results.  This  fluctuation  is  mostly  believed  to  be  the  aftermath  of  the  random  phenomenon 
of  the  seaway  time  history  record  which  had  to  be  modeled  with  a  random  white  noise 
output  process  of  a  digital  filter.  As  a  measure  to  compensate  for  the  fluctuation  and 
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obtain  the  best  results,  the  output  data  were  calculated  fifty  iterations  and  averaged  to  get 
the  final  data  value.  The  graph  in  figure  39  depicts  the  mean  cross-correlation  coefficient 
flanked  by  the  standard  deviation. 

p{t  +  t)  and  p{t  +  t)  were  compared  over  the  entire  time  record,  while  the  following 
expression  was  used  to  formulate  the  cross-correlation  coefficient. 

N 

p(t  +  r)p{t,  t  +  r)l  N 
std{p)std{p) 
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Figure  26.  Forward  Prediction  Diagram 
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Figure  27.  Time  Record  of  P-M  Target  Spectrum  (Y)  versus  Measured  Presure  (p), 

(Predictor). 
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Figure  28.  Zoomed  Plot  of  Time  Record  of  P-M  Target  Spectrum  fy)  versus 
Measured  Pressure  (p),  (Predictor). 
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Figure  29.  Time  Record  of  Measured  Pressure  (p)  versus  Measured  Pressure  p(t+T) 
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Figure  33 .  Time  Record  of  P-M  Target  Spectrum  (T^  versus  Measured  Pressure  (p), 

(Predictor). 
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Figure  34.  Zoomed  Plot  of  Time  Record  of  P-M  Target  Spectrum  (Y)  versus 
Measured  Pressure  (p),  (Predictor). 
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Figure  37.  Time  Record  of  Measured  Pressure  p(t+  t)  versus  Estimate  of  Forward 

Pressure  {t  +  r)) . 


2.4  22.6  22.8  23  23.2  23.4  23.6  23. 

Time(seconds) 


Zoomed  Plot  of  Time  Record  of  Measured  Pressure  pft+rj  versus 
Estimate  of  Forward  Pressure  p(t,  (/  +  t))  . 
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Figure  39.  Mean  Cross-Correlation  Coefficient  and  Standard  Deviation  between  P-M 

Target  Spectrum  and  Digital  Filter. 


V.  CONCLUSION  AND  RECOMMENDATIONS 


A.  CONCLUSION 

In  previous  studies  of  hydrodynamic  effects  on  underwater  structures,  attempts  to 
simulate  waves  using  the  Pierson  Moskowitz  spectrum  have  long  been  of  interest.  This 
work  has  revealed  that  using  an  eighth-order  filter  is  a  convenient  means  of  matching  the 
P-M  spectrum  and  predicting  future  seaway  surface  elevation  behavior. 

It  was  discovered  that  the  application  of  normalized  coefficients  to  model  seaway 
for  various  states  proved  to  be  inconsistent  with  the  actual  behavior  of  sea  waves.  This 
discovery  was  explored  when  a  comparison  of  spectral  densities  and  time  history  plots 
between  dimensionless  and  an  eighth-order  transfer  function  coefficients  was  made.  From 
comparing  the  spectral  densities  and  time  history  of  Spanos’  model  and  the  eighth-order 
filter  with  the  spectral  density  of  the  P-M  spectrum,  it  became  apparent  that  Spanos  model 
failed  to  match  the  P-M  spectrum  while  the  eighth-order  system  proved  to  be  promising  , 
matching  the  P-M  spectrum  very  closely. 

It  has  been  concluded  that  an  eighth-order  transfer  function  provides  good 
coefficients  for  tracking  the  P-M  target  spectrum  and  predicting  seaway  elevation 
response.  What  cannot  be  concluded  fi-om  this  particular  work  is  the  degree  of  optimality 
to  which  the  eighth-order  system  can  actually  match  the  target  spectrum.  Thus  far,  based 
on  history  of  previous  studies  in  this  area,  the  work  conducted  herein  has  indicated  a 
pronvising  direction  for  better  predicting  future  seaway  response.  From  this  research,  it  is 
indicated  that  the  eighth-order  system  allows  adaptivity  and  a  degree  of  future  prediction 
capability  based  on  past  and  current  data. 
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B.  RECOMMENDATIONS 


This  filtering  system  is  merely  a  baseline  version  for  modeling  seaway,  the  gain 
matrix  for  the  ixmovator  was  not  optimized.  That  is,  the  Linear  Quadratic  Regulator  was 
not  optimally  chosen  -  the  weighting  matrices  values  were  based  on  intuition  drawn  form 
past  experiences  without  access  to  experimental  data.  If  the  Q  and  R  matrices  are  chosen 
the  predictor  would  most  likely  predict  further  into  the  future  with  more  accuracy.  Also, 
the  ability  of  the  prediction  model  to  be  adapted  to  seaway  spectra  that  vary  in  peak 
frequency  and  frequency  spread  should  be  explored. 

Also,  it  should  noted  that  the  P-M  target  spectrum  is  recognized  as  an  asymptotic 
form  and  it  has  some  limitations,  for  shallow  water  applications.  Therefore,  the 
Bretschneider  spectrum,  which  belongs  to  a  two  parameter  spectrum  family  that  permits 
period  and  wave  height  to  be  assigned  separately,  should  be  considered  for  fiiture 
research. 
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APPENDIX:  EIGHTH-ORDER  FILTER 


function  [rho]  =  predict (delta) 

Matlab  script  to  function  as  a  Eighth-Order  Digital  Filter  for 
Predicting  Future  Seaway  Elevation  Response 
Tracks  and  matches  Pierson  Moskowitz/ Pressure  Profile 
Spectrum  and  predicts  responses  one  full  wave  length 


Variables 

h  =  Fignificant  wave  height 
w  =  Frequency 

S  =  Defines  Pierson  Moskowitz  expression  for  fully  developed  seas 

t  =  Time  vector 

zeta  =  Damping  ratio 

T  =  Period  interval 

A,B  =  Continuous  plant  model 

Phi, Gamma  =  Discrete  plant  model 

K  =  Filter  Gains 

error= 

x2  =  Estimate  of  state  vector 
Y  =  System^ s  output 

rho  =  Cross-correlation  coefficient 

delta  is  the  number  of  time  steps,  (i.e.  Phi  is  delta  minus  one) 


h=3; 

w=[0.3:0.05:3] ;dw=0.05; 

[1  ,ms ] =size (w) ; 

S=w; 

S  (i)  =8  .  le-3*32 . 2''2  /  (w  (  i )  "5 )  *  exp  (-33 . 56/h^2/  (w(i)  ^4)  )  ; 
end; 

ws= [0-3:0.05:3] ; 


lainbda= [171. 61, 146. 52, 127. 62, 112. 86, 100. 98, 91. 21, 83. 01, 76. 02, 70, 64. 73  . 

60  08  55  94,52.23,48.88,45.83,43.05,40.49,38.14,35.96,33.93,32.05,30.29,... 

2  8. 64  27. 1,25. 6 6, 24. 3, 23. 03, 2 1.84, 20. 71, 19. 66, 18. 67, 17. 75, 16. 88, 16. 07, 15. 30 - 

14.59,13.92,13.29,12.7,12.15,11.63,11.15,10.69,10.26,9.86,9.47,9.11,8.77,  . . . 
8.45,8.15,7.86,7.59,7.33,7.08,6.85] ; 


%  Pressue  Measurement  Simulation 


lambda=  lambda  .'*3.28; 

t=[0:0. 1:200]  ; 

Y=zeros ( 1 , length { t ) ) ; 
for  i=l:length(ws); 

phi=:rand{l,  length (ws)  )  ; phi  =phi -mean (phi)  ; 
y  (  i  ,  :  )  =  (cosh  (  6 *pi  /  lambda  ( i )  )  /  cosh  (40 *pi/ lambda  ( i )  )  )  .  .  . 

* (sqrt {S(i) *2*dw) ) *cos (ws (i) *t+phi (i)  *pi*2 )  ; 

end; 

for  j=l : length (ws) ; 

Y=Y+y { j . : ) ; 

end; 

%  Defines  Coefficients  of  transfer  function  of  Eighth-Order  Filter 


T=0.1;w0=1.28;zeta=0.4; 

n=[l/w0,0] ;d=[ (l/w0)"‘2,2*zeta/w0,l] ;  %  Defines  Fourth-Order  Expression 

num=conv(n,n);num4=conv(num,num); 

den=conv(d,d) ;den4=conv(den,den)  ;  %  Defines  transfer  function 

%  Defines  innovater  gains  for  subject  filter 
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[A1  ,B1,  C1,D1]  =tf2ss  (nuin4  ,  den4  )  ; 

[Phi,  Gainma]  =c2d(Al ,  B1 ,  T)  ; 

K=dlqr  (Phi '  ,  Cl  ^  ,  eye  ( 8 , 8 )  ^*^10 ,  .  1 )  ;  eigl=abs  (eig { Phi-K' *C1 )  ) 
x2=:zeros (8, length  (t) ) ;x4=x2 ;P2=zeros (1, length (t) ) ; P6=P2  ; 

%  For  loop  for  closed  loop  filter 

for  i=l : length (t) 

x2 ( : , i  +  1 )  =  { Phi ) *x2 { : , i) +K' * { Y ( i ) -Cl*x2 { :  ,  i )  )  ; 
x60 ( : , i  +  1)  =  (Phi^delta) *x2 ( : , i  +  1)  ; 

P2 (i) =Cl*x2 ( : , i) ; 

P60(i)=Cl*x60( : ,i) ; 
end; 

for  i=l:(length(t)~delta); 

error60 (i) = (P60 (i) -Y(i+delta) ) ; Y60 ( i ) =Y ( i+delta) ; 
end; 

%  Calculates  cross-correlation  coefficient  between  P-M  spectruia  and  filte 

r 

rho= (Y60*P60 (1:  (2001-delta) )  ' / (std ( Y60 ) *std (P60 ) ) ) / (2001-delta) ; 
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